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HEAT CONDUCTION AS AN EQUIVALENT EIGENVALUE PROBLEM 

by John H. Lynch 
Lewis Research Center 

SUMMARY 

The nonhomogeneous heat- conduction equation has been examined subsequent to its 
transformation into a discrete spectrum eigenvalue form. Theorems have been derived 
which relate the average source temperature to variations in the thermal conductivity 
and/or source over any part of the system. These theorems are applicable to purely 
conducting systems composed of simply or multiply connected sources and nonsources 
having mathematical interface boundary conditions (continuity of temperature and heat 
flux) . It is expected that these theorems will find practical application for conduction 
systems, much like the perturbation formulae of nuclear reactor physics. The theorems 
are general, however, and are not restricted in application to heat conduction systems. 
The application of these methods to any system behavior that can be described by the 
heat -conduction equation requires only a physical interpretation of the eigenvalue. 


INTRODUCTION 

Frequently, in parametric survey and design feasibility studies, one is interested in 
some parameter that characterizes the overall performance rather than the detailed be- 
havior of the system in question. Examples of such parameters would be the effective 
multiplication constant for nuclear reactors or the total heat-transfer rate for heat ex- 
changers. This report describes the development and application of theorems which may 
be used to relate changes in the parameters that characterize an exclusively steady-state 
heat- conduction system (geometry, conductivities, and source strengths) to the average 
source temperature in all source regions. Taking this source temperature as the char- 
acteristic index of the system performance, several areas of practical application are 
possible. 


The theorems derived herein are used most effectively in analyzing the effects of 
parameter variations in systems that have complicated geometries. Depending on the 
geometrical complexity, both analytical and numerical methods have been widely used to 
study the behavior of heat conduction in solids. For simple geometrical configurations, 
the geometry is often approximated by a sphere, cylinder, or slab, and a direct analyti- 
cal solution is obtained. For geometrically complex systems, however, the relaxed 
physical interpretation resulting from approximating the actual geometry as a sphere, 
cylinder, or slab could lead to significant error. Consequently, for complex geometries, 
recourse is usually made to some numerical method. These numerical methods require 
representation of an actual continuous system by the finite difference analog of the con- 
duction equation in discrete space. The finite difference equations are usually solved by 
successive over relaxation or some similar iterative technique. The result obtained from 
such a numerical calculation is unique to the particular geometry, material, and source 
distribution assumed. In order to assess the effects of varying some parameter (mate- 
rial, geometry, etc.), a separate computer calculation is required for each variation. 
Such studies very often require a considerable amount of computer time. The theorems 
developed herein require only a single reference calculation of the temperature distribu- 
tion and appropriate volume integrals (which involve no iteration). These theorems thus 
permit certain types of parameter survey studies to be performed by hand; hence, a con- 
siderable savings in computer time is realized. 

It should be noted that, although the theorems as presented herein are shown to apply 
to heat -conduction systems, this is in no way a restriction. They are mathematically 
general and may be applied to any system behavior that can be described by the heat- 
conduction equation. Their derivation was motivated by a desire to explore the physical 
significance of the eigenvalue of a Sturm- Liouville equation which has a solution that is 
identical to the heat-conduction equation. Since the Sturm- Liouville equation involves 
both the eigenvalue and the eigenfunction (instead of just a function), it is possible to ob- 
tain additional information about a system when both equations have identical solutions. 
This additional information is obtained by performing classical operations on the eigen- 
value equation and by physically interpreting the resulting expressions. 


FORMULATION OF THE EIGENVALUE PROBLEM 

Consider a bounded heat -conduction system D composed of fixed heat generating 
(source) regions and heat dissipating regions (in which sources are vanishingly small). 
The system may or may not be at steady state. Assuming linearity of the conduction 
equation, the temperature at each point in this system may be described by 
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( 1 ) 


V-K(r)vt(r , 0) + q (r ) = p(r)C (r) d* . ^ 0 ) 

P 00 


where 


p density 

Cp specific heat 

K thermal conductivity 

r argument denoting position 

V gradient operator 

t temperature 

9 argument denoting time 

T T » 

q volumetric source strength 

a 

The temperature distribution defined by equation (1) at any particular time 0 may also 
be envisioned as a steady- state eigenfunction solution of the suitably formed three- 
dimensional (space only) Sturm- Liouville equation, which is 

-V • K(r ) Vt(r, 0) = X(0) T(r, 0)t(r, 0) (2) 

where the complementary function, T(r,0), is appropriately chosen. The only restric- 
tion placed on equation (2) is that at least one eigenfunction must be identical to the tem- 
perature distribution obtained from equation (1) for any selected value of 0. The feasi- 
bility of such a formulation is guaranteed by the flexibility allowed in specifying Y(r, 0). 
Integration of equation (2) over the system volume gives 


- f V-K(r)vt(r,0)dV= X(0) / T(r, 0)t(r, 0)dV 

‘'tn *'T\ 


( 3 ) 


Since the heat flux K(r)Vt(F, 0) is continuous for all 0 (except possibly zero), the left 
side of equation (3) may be transformed by Gauss's theorem to surface integrals of the 
heat flux over all parts of the system and is recognized as the total heat loss rate from 
all parts the system at time 0. Thus, 
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V • K(r ) Vt(r, 0)dV = - /K(r ) Vt(r, 0) • dS 

1 ^ i 

all | all i 


= system heat loss rate at time 9 


= Me) f T(r, 0 )t(r, 0)dV 


( 4 ) 


where £ denotes some part of the system, and the sum of the £ equals D. Normally, 
only the boundaries of the system contribute to the divergence integrals. The summation 

A 

has been used to allow for possible discontinuities in the heat flux at 0 = 0. The eigen- 
value is now defined as the ratio of the sum of the heat loss rates to heat production rate 
for the entire system at time 0: 


yW)Vt(r,0)-dS 


\(0) = _allj_ 


f q'(r)dV 
u 


Losses/sec 

Productions/sec 


This definition fixes the integral of T(r, 0)t(r, 0) as 


( 5 ) 


/ T(r, 0)t(r, 0)dV = f q (F)dV = total system heat source (6) 

D u 

but places no further constraint on T(r, 0), so that the formulation of equation (2) is still 

A X 

possible. The form of T(r, 0) is not known for all values of 0, however, it can be seen 
that, for the steady-state (0 = «), A is 1. 0 and the functional form of T(F, «) must be 

T(r,°^ = q(r)/t(r,«^ (7) 

This guarentees the equivalence of equations (1) and (2) for the steady state but is not 
generally valid at all values of 0 since each point may be out of balance by a different 
amount. The steady-state form of this complementary function, given by equation (7), 
thus transforms the nonhomogeneous steady-state conduction equation into a homogeneous 
eigenvalue form with known linear operators. 
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VIRTUAL IMBALANCE IN THE EIGEN SYSTEM 


The concept of X, introduced in the preceding section as a measure of the time- 
dependent ratio of net losses (Btu/hr) to productions (also Btu/hr), permits a physical in- 
terpretation of the heat equilibrium of the system. If X is greater than 1. 0, the system 
is cooling down, and if X is less than 1. 0, it is heating up. For X equal to 1. 0, the 
system is at equilibrium. The definition of X as an eigenvalue also enables the develop- 
ment of a method for determining the degree of virtual imbalance for the equivalent 
eigensystem when the physical properties have experienced some alteration. The virtual 
imbalance referred to herein is the change in the eigenvalue of the equivalent eigenvalue 
equation (eq. (2)) when the operators are perturbed. It is not a function of time since the 
eigenvalue equation is not time dependent but is associated with a particular time. Ex- 
amples of such alterations would be a change in material at a particular location, the in- 
troduction of additional sources into the system, or changes in dimensions of subregions 
of the system. These changes may all be effected by changing the spatial distribution of 
the thermal conductivity and/or source functions. The following derivation is concerned 
with the virtual imbalance in the eigensystem, designated as 6X, associated with alter- 
ations in the conductivity in any region of a geometrically arbitrary system. The ap- 
proach used is similar to the derivation of first-order eigenvalue changes by the theory 
of small perturbations. It is shown in detail to emphasize the inherent differences from 
perturbation theory that arise from the requirement that the eigenfunction must also be 
the solution to the nonhomogeneous conduction equation. 

For simplicity, the system is assumed to be at the steady state. For this condition, 
the 9 argument is dropped, and equation (1) becomes 

-V • K(r ) Vt(r) = q'(r) (8) 

For this case, the eigenvalue of the eigentemperature distribution defined by equations 
(5) and (7) is 1. Q-and is denoted by Xq. Equation (2) is then 

-V • K(r) Vt Q (r) = X 0 T(r)t Q (r) (9) 

To facilitate the manipulations that follow, this is rewritten in operator form 

IT 0 = V T 0 <10 > 


where 


I s -V • K(r)v 


( 11 ) 
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J = T(r ) 


( 12 ) 


Tq h (zero subscript denotes the unperturbed steady state) (13) 


Next, a function is considered that will be called the adjoint eigentemperature. In 

jjc 

the discussion that follows, this will be designated by Tq. This function is defined by 


* * — * * 

1 T 0 - V T o 


(14) 


and has the same boundary conditions as the solution of equation (10). 

If the system experiences a change in thermal conductivity in some subregion R 
contained in D, the loss operator I of region R will change by some amount 61. The 
eigenvalue will change by some amount 6A and the eigentemperature distribution by 6T. 
In order that the source remain unchanged, J must be altered by an amount 6J such 
that 



(J + 6J)(Tq + 6T)dV = 




T t t 

q (r)dV 


(15) 


The defining equation for the perturbed system is then 

(I + 61) (Tq + 6T) = (Aq + 6A)(J + 6J)(Tq + 6T) (16) 

Multiplication of the perturbed equation (eq. (16)) on the left by the unperturbed adjoint 
$ 

eigentemperature Tq, and the unperturbed adjoint equation (eq. (14)) by the perturbed 
temperature (Tq + 6T) gives 

Tj(I + 61) (Tq + 6T) = (Aq + 6A)Tq(J + 6J)(Tq + 6T) (17) 


and 


tJi*(Tq + 6T) = Aq(Tq + 6T)J*Tq (18) 

Subtracting equation (18) from equation (17), integrating the result over the system vol- 
ume, and making use of the identity, 
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(19) 


f<PS dV=/ g 


g(p * f dV 


(where f and g are any real functions satisfying the same homogeneous boundary con- 
ditions, and <p is a linear operator) gives, for 6A, exactly 


f T*6I(T 0 + 6T)dV - (X Q + 6A) f T*6J(T Q + 6T)dV 

6A = _R Q 



+ 6T)dV 


( 20 ) 


where Q denotes the source region. A significant point here and one worthy of special 
emphasis is the exact property of equation (20). In the conventional perturbation theory 
development of eigenvalue perturbations, recourse is made to ignoring products of differ- 
ential magnitude thus obtaining an approximation for 5 A and restricting application of 
the theory to small perturbations exclusively. No such procedure, however, is required 
to obtain equation (20); thus, the range of applicability is unlimited for certain types of 
problems, as will be seen. There are several approximate forms of equation (20) that 
may be obtained by using T Q in place of (T Q + 6T) in one or more of the integrals. If, 
however, in the interest of preserving the exact nature of the development, this limit is 
taken in all integrals, the second integral in the numerator vanishes by equations (15) 
and (6). The resulting expression is 


/ T O 6IT 0 dV 

6A = -^ = lim (6A) (21) 

r * 6T-0 

T 0 JT 0 dV 

This may be evaluated in general (for zero or nonzero boundaries) using equation (58). 
The admissability of the limiting operation is a point that may be arbitrary for the ana- 
lytic version of equation (21); this was not explored. However, for numerical represen- 
tations of the integrals (the primary source of application of the theorems to be derived 
is numerical) in equation (21), this operation is correct. It should be noted here that 
this limit is not the same as the limit that would be obtained using the definition of A 
given by equation (5) at the instant of the perturbation In taking the limit, another de- 
viation from perturbation theory had been introduced. The definition of A is the same 
initially (before the perturbation) whether the weighting function is present or not as can 
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be seen by comparing equations (5) and (24). However, the limits of 6A as 6T ap- 
proaches zero are not. Equation (21) then describes the virtual imbalance in the eigen- 
system that would result from a change in conductivity in some or all parts of the sys- 
tem. 

It is of interest to derive equation (21) by using variational formulation in order to 
demonstrate that identical results may be obtained with an alternate approach. In oper- 
ator form, 


IT o = V T o 


(22) 


and 


jfc jje jjc 

I Tq = V Tq (same eigenvalue) 


(23) 


so that a variational principle (Rayleigh quotient) for Aq is (ref. 1) 


X 0 = 


^ T o IT o dv 


dV 


^V T o dv .^vX 


dV 


(24) 


Varying equation (24) gives, after trivial manipulation, 


/ T o 5IT o dv - x o / T 0 MT 0 

5A = -5 5 


dV 




Q 


T o JT o dv 


(25) 


Again, if the limit is taken as 6T approaches zero, no variation in J is required. 
Thus, 6J is zero and equation (25) becomes 


6A = 


/ T 0 6IT 0 


dV 




(26) 


T 0 JT 0 dV 
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which is identical to equation (21). The perturbed eigentemperature (Tq + 6T) is of no 
interest since it has no known physical significance. This arises from the fact that J is 
considered linear (independent of T) throughout the perturbation and loses contact with 
the physical problem for any temperature other than the one for which it was specifically 
chosen. It should be observed that, if the perturbation is made over most of Q(i. e. , if 
R - Q), 6 A may be approximated by relatively crude estimates of the temperature dis- 
tribution, the error in the 6 A so calculated being much less than that of the tempera- 
tures used for evaluation (ref. 1). This is because when R =* Q, equation (26) is a sta- 
tionary form for 6 A; that is, 6 A in this case may be taken as the eigenvalue of the 
equation 


6IT q = 6AJT 0 (27) 

the Rayleigh quotient of which is formed by pre multiplication by Tq and integration over 
the volume giving equation (26). This may afford some computational convenience when 
evaluating equation (26) for system changes of this type. 

It should be noted that the limiting form 6 A does not involve J or 6J, so that the 
functional form of T is not required to evaluate 6A. Hence, the perturbation can be 
made in a system that was not necessarily at the steady state initially. This was as- 
sumed for clarity of presentation but is not a requirement. Hence, the zero subscript in 
all of the preceding derivation may refer to an initial distribution that has not reached 
the steady state. 


TEMPERATURE PREDICTION - THEOREMS I AND II 

A utilitarian aspect of the eigenvalue approach to heat conduction is the prediction of 
the change in average source temperature resulting from a change in conductivity and/or 
source anywhere in the system. In the following discussion, the J coefficient operator 
is not fixed throughout the perturbation derivation as was done previously. Instead, the 
time -dependent nature of X is considered. In the previous discussion, the operator J 
was taken to be fixed and associated with a particular temperature distribution and time; 
thus the perturbed eigentemperature was meaningless. For the purposes of the develop- 
ment that follows, this operator is thought of as varying continuously at all times follow- 
ing the perturbation in such a manner as to ensure that the real temperature and eigen- 
temperature remain identical. 
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Conductivity Effects 


Consider the steady-state form of Aq described by equation (24). The ratio Aq is 
understood to be an eigenvalue in this equation. If the time dependence of the tempera- 
ture and the J operator is recognized, the eigenvalue variation will also be time de- 
pendent. The eigenvalue variation is given in operator notation by equation (20), which, 
for the purpose of this derivation, will be considered to incorporate a time dependence in 
T, J, and A. Following a conductivity perturbation, 6A will assume some transient 
functional behavior. In the limit as 9 approaches infinity, 6A must become zero since 
the system tends to regain the steady-state condition (X approaches 1.0). This means 
that 6J ultimately takes on some nontrivial form in view of the requirement of a con- 
tinuously varying J operator just discussed. Denoting the final temperature distribution 
taken on by T + 6T as T , the perturbed steady-state distribution (eq. (20)), becomes, 

— ~ Jr 

in the limit as 6\ goes to zero, 

/ T o SIT P dv - < x 0 - 5X) / T o 6 JT p dv 
5A = -5 Q = o (28) 

dv 


so that 


/ 


V IT P dv 


\> T J 6jT P ™ 


J T A JT„ dV 


* 

vs 


/ T a JT„ dv 

•'n 


* 

V"p 


(29) 


The adjoint and real eigentemperatures are identical, since equation (10) is self adjoint 
(ref. 2). The self adjointness of the homogeneous Sturm- Liouville equation only holds 
for homogeneous boundary conditions of the form cxd Tq/9F + 0Tq = 0. For the present, 
these conditions are assumed, that is, zero temperature, zero heat flux, or zero exter- 
nal sink temperature. The treatment of nonhomogeneous boundary conditions is left for 
subsequent discussion. Since the system is initially at steady state, Aq is 1. 0. Re- 
placing Tp by Tq in the denominator only of equation (29) (the denominator is arbi- 
trary) and dropping the adjoint distinction give 
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( 30 ) 


/ T n 6IT dV / T n 5JT_ dV 

‘'d v P */r\ 


_ Q 


0 UUi p 


^ T o JT o dv 4 VT ° 


dV 


Multiplying the integrands of equation (15) by Tq gives 

T 0 ’q = T q (J + 6 J)(T q + 8T) = T Q (J + 6J)T p 


(31) 


or 


T 0 6JT p - To? - T 0 JT p < 32 > 

From equations (8), (9), (12), and (13) with Xq equal to 1.0, 

T 0 JT q = qT* (33) 

Substituting equation (33) into equation (32) gives 

T 0 5 JT p= q( T o- T p) (34) 

The 61 operator will be seen in a subsequent section to involve only the temperature 
gradient. Hence, T p may be replaced by Tq exactly in the numerator of the left side of 
equation (30) if the integrand is multiplied by the ratio of the unperturbed to perturbed 
conductivities in region R. This is only valid, however, when the heat flow to all points 
in the perturbed region remains constant (from Fourier's law, the temperature gradient 
at a point is inversely proportional to the conductivity for constant heat input). If the 
heat source to the perturbed points changes as a result of the perturbation, a different 
correction may be required. This is discussed in appendix A. The former correction 
will be retained herein to simplify the derivation. 

Making this replacement and substituting equations (33) and (34) into equation (30) 

give 



J Q '< i'< T o - V dv 


f q’ T o dv 

J Q 0 


(35) 


11 



Comparing equations (35) and (21) it is noted that, if the correction K ^/Kq is taken out- 
side of the integral as would be possible if the conductivities were constants over R, 
equation (35) could be written 



f »" 
J * 

Q 


(Tq - T > dv 


A 

J Q 


T o dv 


(36) 


where 6A is given by equation (26). Also, if the source strength is spatially constant, 
the numerator and denominator of equation (36) may be divided by the source volume and 
rearrange to obtain the fractional change in the average temperature of the source due to 
a conductivity perturbation anywhere in the system. 

Theorem I 



This theorem relates the fractional temperature change in the source region to the 
limiting form of the eigenvalue change that is caused by an alteration in the loss operator 
of the equivalent eigensystem. The conductivity ratio has been assumed constant over R 
in this derivation so that equation (37) could be shown simply. Making this assumption is 
not necessary, however, because, for any numerical model, conductivities must be 
treated numerically as discontinous step functions. Even in practice, most systems are 
considered to be composites of materials each having its own conductivity. If the con- 
ductivity is temperature dependent, equation (1) is nonlinear; this case is not amenable 
to the eigenvalue approach, since the operators must be assumed linear. If the source 
is not spatially constant, a different type of average change is obtained, namely, the 
change in the source weighted temperatures. The limiting 6 A is evaluated using 
equation (57). 

An interesting corollary implied by equation (37) is that the relative change in the 
source temperature resulting from two different changes in the conductivity in a specific 
region is inversely related to the conductivity ratio and directly related to the ratio of 
the conductivity changes: 
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Corollary to Theorem I 



( 38 ) 


This is derived by dividing equation (37) for one change (subscript 1) by equation (37) 
for a second change (subscript 2) in the same region. Then, by equation (57) the tem- 
perature distribution dependence may be removed by canceling. 

Equation (38) enables comparisons of the relative effects on source temperature of 
using several different materials in a specific region of the system with no knowledge of 
the unperturbed system geometry, materials, or source strength (only the unperturbed 
conductivity in the region of interest is required). 

In all of the foregoing, the initial state of the system was assumed to be steady. 

This is not an essential assumption. The Aq could have been carried throughout the de- 
velopment thus giving the fractional change in source temperature based on the average 
source temperature at the time of the perturbation (which could have been changing with 
time) and the steady-state temperature finally achieved. This has little practical signi- 
ficance, however, and will not be given detailed consideration. 


Source Effects 


To derive an expression for the effect of changes in the source distribution on the 
unperturbed average source temperature, manipulations shown in equations (28) to (36) 
are performed. This gives, when it is remembered that the integration is over the 
whole system instead of with the limits shown in equation (28), 


6 \ = 



T*5J(T Q + 6T)dV = 0 


(39) 


This expression reflects both the change in source distribution because of the perturba- 
tion and the change in J in the unperturbed (original) source that is required to offset 
the temperature change that occurs as a result of the perturbation. Distinguishing these 
regions by their integral limits, equation (39) may be written as 



T 0 6J(T q + 6T)dV = - 



T q 6J(T 0 + 6T)dV 


(40) 
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or 



dV 


-/ 


Ml 

q(T - 


T n )dV 


(41) 


The last step is similar to manipulations performed in the previous derivation. Reten- 
tion of the minus sign depends on how 6q is defined. On the left of equation (40) , use 
has been made of the fact that, prior to the perturbation, J was related to q by 

JT q = q (42) 

If the perturbed region contained no source originally, J was zero and, after the pertur- 
bation, 


6J(T q + 6T) = 6q' (43) 

If the sources are spatially constant, manipulation of equation (41) and replacing Tq 
with t^ give 

Theorem n 


6T _ 6q 

nn in 

q 


dv 


(44) 


The 6q is shown outside of the integral with the numerical model in mind. If the 
change were not spatially constant, it would have to be represented by a set of constant 
values over small regions. The total effect would then be found by applying this equation 
to each small region. This is also true for theorem I. Theorem II relates the fractional 
change in source temperature (the source being defined as all regions containing sources 
before the perturbation) to the addition or deletion of a source in any region R. 

The corollary to theorem II is formed as was the corollary to theorem I. This is 

Corollary to theorem II 

6T\ (fiqV 

= -111 (45) 

6T 9 Ia'Z\ 

2 (Sq^ 
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As with the corollary to theorem I, this result is not obvious. It allows a direct compar- 
ison of the average temperature increases in the source region of a system, which would 
result from the addition of sources of varying strength in a given location, with no know- 
ledge of the system properties or any temperature redistribution that may result from 
the introduction of an extraneous source (only the unperturbed source in the region of in- 
terest is required) . 


EVALUATION OF 5A FOR CONDUCTIVITY OR SOURCE CHANGES 

In order to emphasize some of the physical restrictions that are encountered in the 
application of equation (21), its algebraic equivalent must be shown in detail. 

The denominator of equation (21) presents no problem since, by equations (8), (12), 
and (13), it is just 


Denominator = 



Mt — . 

t Q (r)q (r )dV 


(46) 


The integrand of the numerator is derived as follows: 

1= -V-KV 

I + 61 = -V- (K + 6K)V 


(47) 

(48) 


By identity, 


V • SV = VS • V + SV • V (49) 

where S is any scalar and V is any vector. Using equation (49), equation (48) be- 
comes 


1 + 61 = 



+ 6K) • V + (K + 6K)V 


2 ] 


(50) 


Also, equation (47) can be written as 

1= - (VK- V + KV 2 ) (51) 

Subtracting equation (51) from equation (50) gives 
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61 = - (v5K • V + 6KV 2 ) 


(52) 


Substituting equations (52) and (46) into equation (21) and changing T Q and T Q to their 
equivalent t^ give 


6A = 


- ^(t 0 * 5 K-*o + t 0 5K A>) dV 

t t T 

qt Q dV 



(53) 


Making use of the same identity as before yields 

V • (t 0 6KVt 0 ) = Vt Q 6K • Vt Q + t Q 6KV 2 t 0 

= 6K | Vt Q | 2 + t 0 6KV 2 t Q + t 0 V6K • Vt Q 
Solving for -(tgV6K' Vtp) results in 

-(t Q v6K • Vt Q ) = -V • (t Q 6KVt 0 ) + 6K | Vt Q | 2 + t Q 6KV 2 t 0 
Substituting equation (54) into equation (53) gives 


6A = 



(54) 


(55) 


At this point it is recalled that the boundaries of the system were assumed to be at 
zero temperature. Then, by Gauss's theorem, the first integral in the numerator of 
equation (55) can be transformed to a surface integral over the surface S of D to yield 

- [v • (t 0 6KVt 0 )] dV = f - [v • (t 0 6KVt 0 )]dV = - ./t 0 6KVt 0 • dS = 0 (56) 

Hence, for the zero surface temperature condition, 
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( 57 ) 



If the boundary temperatures are nonzero, the application of equation (57) in equation (37) 
requires that the temperature be some constant at the physical boundaries of the system. 
It can be seen that, for this case, the gradient is unaltered and the denominator is in- 
creased by the constant surface temperature. The fractional temperature change in the 
source 6T/T is affected identically by the addition of a constant to the temperatures in- 
volved. Hence, this is only a normalization effect and equation (57) remains effectively 
valid for use in equation (37) to predict fractional average source temperature changes. 

The 6A computed using equation (57) does not represent the virtual imbalance of an 
eigensystem with nonzero boundaries. If the boundary temperature is some constant, 
however, t^ in equation (57) is taken to be the solution to the homogeneous boundary 
value problem formed by subtracting the boundary value from the temperature (ref. 3). 
All operations required to derive equation (57) may be performed exactly as shown on 
this temperature excess function (which has been normalized to zero at the boundaries) , 
and equation (57) will again be obtained. Hence, to use equation (57) for the nonhomoge- 
neous boundary condition problem (prescribed nonzero constant), it is noted that the ad- 
dition or subtraction of a constant to the numerator results in no change. However, the 
denominator must be decreased appropriately. The exact expression for the nonzero 
boundary case is then 


6A = 



(58) 


This correction is never used, however, since it is not required for the 6A that appears 
in equation (37). This is because the boundary temperature must also be included in the 
denominator of the change in the average source temperature. This is apparent in the 
simple example calculations that follow. The 6T/T is readily seen to be independent of 
the boundary temperature by adding a constant to the temperature in the conventional so- 
lution. Equation (44) does, however, require slight modification when used with a con- 
stant nonzero boundary problem. For this case, the numerator should be corrected to 
give 
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( 59 ) 


6T _ R 
T 


J 5qt 0 dV - US) J 6q dV 
R u U ^R 


f q t 0 dv 

y Q u 


EXAMPLE CALCULATIONS 


In order to demonstrate the mechanics of the temperature prediction technique 
(theorem I) proposed herein, a sample calculation is performed for an elementary geom- 
etry. The average temperature change in the source region for a simple one-dimensional 
system is calculated for conductivity perturbations in source and nonsource regions. 
Theorem II is demonstrated by using a source perturbation in the same geometry. These 
calculations are elementary examples of the application of the 6A method and serve the 
dual purpose of illustrating the method in detail and providing test cases necessary to 

complete the exposition. They are trivial examples 
in the sense that the results are obvious and may be 
obtained more directly. They are not suggested as 
typical applications, but are included only to help 
clarify the concept. To examine theorem I, the infi- 
nite slab at steady state (see sketch) is considered. 

The conductivity in region 1 is to be perturbed. First, 6T/T will be computed 
using theorem I. Denoting regions by subscripts and perturbed and unperturbed states 
by superscripts, a conductivity perturbation in region 1 yields 


= o 

dx x=0 


t'At 


Non source 


Region \A Re 9 ion 2 


T(L + H} - 
X 



(60) 


The temperature distributions assuming continuity of heat flux and temperature at x = L, 
symmetry at x = 0, and zero at x = L + H are 
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(61) 


,TT Ml 

:? = — 1- (l 2 - x 2 ) + ± 

A n 


2K, 


LH 

k» 


H! T t t 

T° = -9- (l 2 - XL) +-3Jl5 


k; 


k; 


(62) 


The gradient in region 1 is 


ST” 

0X 


t T T 

_ q x. 

rf’ 


ST, 


Sx 


"’2 2 

= q x 

4 


(63) 


which gives 



and the fractional change in the source temperature is, by equation (37), 


/ \ 


T 0 kP 



(65) 


This is compared with the result obtained conventionally as follows: 6T/T in the source 
computed using equation (61) gives 
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( 66 ) 


and 



(67) 


( 68 ) 


which is identical to the result obtained in equation (65). It is noted that, for a conduc- 
tivity perturbation in the nonsource region, both methods give 


6T 

T 



(69) 


20 


* In both of these calculations, the Kq/K^ correction was used since the heat flow to the 
perturbed regions was not affected because of the symmetry of the perturbation. 

To minimize notation, theorem n may be readily examined by assuming that the con- 
ductivities are the same in all regions. Let the perturbation consist of adding a source 
to region 2 having the same unit strength as in region 1 before the perturbation. The 
change in source strength in region 2 is then 


f 1 1 t t T T 1 1 A t t t t » t t T f 

6q= q 2 q 2 =q i~ 0=q i sq 

Theorem n gives, for the fractional change in T^, 


ST 

T 


6 ’q /t 0 dV 
R u 

-4' <1 ’ T o dv 


/’L+H 0 

L 1 2 dx 

/ L T»dx 
Jq 1 


(70) 


(71) 


ST _ 3H 2 

T 2L 2 + 6LH 


(72) 


It t 

After the perturbation, q will be in both regions. Since K was assumed the same in 
both regions, we have a homogeneous slab of thickness L + H. The temperature distri- 
bution in the perturbed case is then 


? ft p 

= (L + H) 2 - x 2 
2K L J 


and the average over region 1 is 


T = -3— 
p 2K 


(L + H) 2 - — 


Using equation (61), the unperturbed average is 


T n = -5L (- L 2 + 2LH| 
0 2K \3 


(73) 


(74) 


(75) 
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giving a fractional change of 


6T _ 3H 2 

T 2L 2 + 6LH 


which agrees with equation (72). 


APPLICATIONS 

It is not the intent of this report to elaborate on the application of the proposed theo- 
rems since the a priori practicality of any concept is difficult to assess depending 
strongly, as it does, on individual ingenuity. It is worthwhile, however, to devote some 
discussion to possible areas of application for the sake of completeness. 


Survey Calculations 

As suggested in the INTRODUCTION, this method is most powerful as a practical 
tool when applied to conduction problems having complex geometries. Suppose, for ex- 
ample, that a design survey is to be performed for some heat -conduction device that, out 
of necessity, has a complicated geometrical configuration. It is also assumed that the 
geometry is of sufficient complexity to preclude any hand calculational approach and that 
the principal interest is in the source temperature. In this situation, the required cal- 
culations must be performed numerically. The disadvantage here is that, as some pa- 
rameter is varied, say the size of a particular part of the system, several machine cal- 
culations must be performed in order to establish a trend. This could be expensive de- 
pending on the complexities of the system and the machine code used to analyze it. The 
eigenvalue method, however, would require only one calculation of the temperature dis- 
tribution for a reference case and an additional calculation of the volume integrals of the 
gradients squared (which would require negligible machine time compared to the tem- 
perature calculation). With this information, several types of parameter variation com- 
putations could be performed by hand using theorem I. 


Nonlinearities 

Effects of nonlinearities are also amenable to hand calculation using theorem I. 
Suppose that, because of extreme temperature variations, the thermal conductivity in 
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the above problem varies strongly with position. To approximate the average source 
temperature, the temperature distribution is first calculated assuming temperature inde- 
pendent conductivities, and then the deviations in the source temperature due to the non- 
linearity in K are summed. 


Equivalence of Perturbations 

Another possible application of theorem I is the calculation of changes in one region 
that are required to offset changes that are made in another region in order for the 
source temperature to remain constant. 


Comparison Calculations 

If the relative effect of several different materials in one specific region is of in- 
terest, the corollary to theorem I may be applied directly. In fact, it requires no de- 
tailed knowledge of the system (geometry, source strength, etc.); only the unperturbed 
conductivity in the region of interest is required. 


Source Interactions 

Consider a large source such as an electrical bus bar embedded in a thick insulator. 
A second source of the same unit strength is to be inserted in the same insulator. The 
effect on the first source (average temperature change) is seen from theorem II to be 
equal to the ratio of the unperturbed temperature integral over the region that would be 
occupied by the second source to this integral over the first source. Hence, if the un- 
perturbed temperature distribution is known, the optimum position (and shape) of the ad- 
ditional source may be selected intelligently so as to minimize the interaction. With 
some ingenuity, several applications of this kind are possible. 


Perturbations 

For problems of complicated geometry, the classical application of perturbation 
theory suggests that a simpler problem be solved which closely resembles the problem 
in question. The solution of the complicated problem is then obtained by perturbing the 
simpler problem. The eigenvalue methods of theorems I and n may also be applied in 
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this manner but have the added advantage of imposing no constraint on the magnitude of 
the perturbation (for problems in which the perturbation does not affect the heat sources) , 
which is an inherent weakness of perturbation theory. 

To summarize the discussion of practicality, it may be said that the practical as- 
pects of these methods are in essence all founded on the ability to anticipate the effect of 
a change in any conduction system using only the unperturbed temperature distribution. 
For geometrically complicated problems, especially those that must be solved numeri- 
cally, the savings in computer time could be significant. 


DISCUSSION OF THE EIGENVALUE METHOD 

It is of interest to note a few of the mathematical properties of the eigenvalue ap- 
proach to heat conduction. Fortunately, the complementary function never requires 
evaluation and acts only as a mathematical expedient. In the derivation of theorem I, in 
effect (Tq + 6T) is replaced by Tq exactly and not approximately as in the first-order 
theory of small perturbations. The significance of this is that theorem I is valid for 
changes of any magnitude for problems in which the heat source to the perturbed points 
is constant, whereas first-order perturbation theory is valid only for small perturba- 
tions. Thus, theorem I is a more powerful calculational device for problems of this type 
than perturbation methods are for eigenvalue problems. For problems having variable 
heat source to the perturbed points, the method can be used as an approximation similar 
to first order perturbation theory. This is discussed in appendix A. 

In addition to this, for R ^ Q, the eigenvalue method described herein retains the 
capacity of yielding an approximation for 6A that will generally be much more accurate 
than the temperature distributions used to evaluate it because of its stationary property. 
Hence, crude estimates of the temperature may be used to obtain relatively accurate 
values of 6A. 


CONCLUDING REMARKS 

The solution of the nonhomogeneous heat -conduction equation for an arbitrary geom- 
etry has been viewed as an eigenfunction over the same geometric space having homoge- 
neous boundary conditions. An approach similar to first-order perturbation theory was 
used to examine changes in the eigenvalue that would be caused by changes in the conduc- 
tivity or source strength distributions. As a result of this, theorems have been derived 
that relate the average source temperature to changes in thermal conductivity or sources 
at any or all points in the thermal conduction system. It is expected that these theorems 


24 


may find practical application for conduction systems much like the perturbation formu- 
lae of nuclear reactor physics. They are especially suited to geometrically complex 
systems that must be analyzed numerically in that they enable certain types of parameter 
variation studies to be performed by hand using the computer solution for a single refer- 
ence case. These theorems were developed out of a recognition of the physical signifi- 
cance of the eigenvalue of an equivalent Sturm- Liouville solution. Their extended appli- 
cation to other physical systems, which can be described by the heat- conduction equa- 
tion, requires only a physical interpretation of 6A. 

Lewis Research Center, 

National Aeronautics and Space Administration, 

Cleveland, Ohio, December 2, 1966, 

122-29-05-08-22. 
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APPENDIX A 


REPLACEMENT OF THE PERTURBED GRADIENT BY THE 
UNPERTURBED GRADIENT IN EQUATION (30) 

Conduction calculations performed numerically yield the temperature at a set of node 
points. In order to obtain the gradient at a point (the gradient is required to evaluate the 
6 A expression), a seven-point-difference description might be used in rectangular co- 
ordinates, for example, where 



j j 




j = x,y,z 


(Al) 


where the unit vectors have been omitted for simplicity and A denotes the appropriate 
mesh interval. Hence, the gradient is seen to be numerically a simple function of the 
temperature differences between i and its immediate neighbors. For simplicity, only 
one of these temperature differences is discussed. A correction to this AT that will 
account for a change in the conductivity of the material separating the two reference 
points is examined. The function required to evaluate the correction might possibly be 
obtained numerically along with the reference calculation of the temperature distribution. 
This is not suggested, however, since a simpler method may be possible depending on 
the nature of the problem and the numerical procedure used to solve it. It is shown 
herein primarily to give some insight into the behavior of the gradient correction and to 
provide qualitative bounds on the range of validity of the temperature prediction method 
when used as an approximation (without correction). 

From elementary, direct-current theory, the voltage drop (which is analogous to 
the temperature difference) across any resistance R is given by 


AV = 



(A2) 


/ / 
where Erp is the Thevenin equivalent electromotive force and Rrp is the Thevenin 

equivalent resistance of the system measured across the points in question in the ab- 
sence of R and any source. 

The ratio of perturbed to unperturbed voltage drops and, thus, temperature differ- 
ences by analogy, across the same points, is then 
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(A3) 


A 

AT p _ % (&? + <%) _ \&q l 

AT 0 (Rq (tflrji + (Rp) / (Rrp \ 

v 1 / 

where the (R values refer to thermal resistances 

<R s (A4) 

j*K 


in which i is the separation distance between the points and s/ is the area of heat 
transfer. 

Qualitatively, a perturbation can be characterized by the volume over which it ap- 
plies and the amount by which the conductivity has been changed. Some insight into the 
behavior of the correction for limiting combinations of these characteristics can be 
obtained using (A3) if the problem is such that the following may be assumed: 


(a) perturbation over a small volume =► (ftp »> (Rq 

(b) small perturbation in K =► (Rq « (R 

(c) perturbation over a large volume =► (Rrp <« (Rg 

(d) large change in K =*► (R p «< <R Q or the reverse. 
Considering all combinations of the above, we have 

Combination 

Large volume, small change in K 
Large volume, large change in K 
Small volume, small change in K 


or (Rp 
or (R p 


A V AT 0 

a 1.0 
a 1.0 
a 1.0 


Small volume, large change in K = y ^ - v k p 

/ 

If the Thevenin equivalent resistances can be computed, an exact correction is pos- 
sible. If these are not available, however, the above considerations should provide 
some feel for the validity of the result. 
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APPENDIX B 


SYMBOLS 


s/ 

area of heat transfer 

r 

vector denoting position 

D 

limit of integration denoting 
whole system 

S 

surface limit of integration de- 
noting surface of whole system 

e t 

Thevenin equivalent electro- 

T 

operator form of t 


motive force 

T 

average source temperature 

f, g 

any real functions 

t 

temperature 

I 

loss operator 

X 

space variable in example cal- 

i 

subscript denoting node point 


culation 


in numerical conduction 
model 

<M 3 
6 

arbitrary constants 
variation operator 

J 

operator form of T 

5\ 

perturbation in A. when X is 

j 

summation index denoting co- 
ordinate axes in rectangu- 


taken to be an eigenvalue 

K 

L, H 

lar coordinates 
thermal conductivity 
region dimension in example 

6A 

X 

limiting of 6X 

ratio of losses to sources for 
whole system 

l 

calculation 

distance between two node 

4 

limit of integration denoting 
some part of system 

points 

T 

complimentary function 

P 

superscript or subscript de- 


arbitrary linear operator 

Q 

noting perturbed state 
limit of integration denoting 

* 

superscript denoting adjoint 
function 


source region 

V 

gradient operator 

t f i 

q 

volumetric source strength 

0 

superscript or subscript de- 

R 

limit of integration denoting 
perturbed region 


noting unperturbed or steady 
state 

Rrp 

(R 

t 

Thevenin equivalent resistance 
resistance to heat flow 

d),(2) 

limits of integration in example 
calculation 
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